#!/usr/bin/perl

# Run in background sequentially:
# (a.out ; b.out ; c.out) &

# Seed the random generator with a unique seed, every time sweeps are prepared
srand (time ^ $$ ^ unpack "%L*", `ps axww | gzip -f`);

# Sweep selective coefficient models
# RSM model 2 is s=0.6 for all 5 loci
$rsm[0] = 4;
$rsm[1] = 5;

open(in, "run.simulation.coefficients.models.csv");
while(<in>){
    chomp;
    @tok = split(/\,/);
    $rsm_loci_string[$tok[0]][$tok[1]] = $tok[2];
}
close(in);

$outputdir = "/home/rumen/simulation1";
if (!(-d $outputdir)) {
    mkdir($outputdir, 0777);
}
chdir($outputdir);

$args = "rm -fr $outputdir/job* $outputdir/infile*";
system($args);

# Grid size
$grid_size_x = 300;
$grid_size_y = 300;

# Sweeps replicates
#$replicates_per_sweep = 1;
$start_replicate = 2;
$stop_replicate = 2;

$division_rate = 0.004; 
$division_max_rate = 1;
$division_min_rate = 0.001;

$death_rate = 0.003; 
#$death_max_rate = 0.1;
#$death_min_rate = 0.005;
#$death_rate = 0.005; 
$death_max_rate = 0.1;
$death_min_rate = 0.00025;

# Days interval to dump grid state, clone frequencies, clone patterns
$dump_grid_state_t = 182;

# Days to cleanup phylogeny, record number of extant clones
$cleanup_t = 1000;

# Days of simulation (20 years)
$simulation_t = 14600;

# Number of microsat/neutral loci
$neutral_loci_n = 100;

# Neutral microsatellite mut rate is 10^-5, 10^-6, 10^-7 (literature based)
# Run with higher neutral mutation rate 1-e4
$nu[0] = 0.0001; # 1e-4

# Mutation rate of the selective loci 10^-6 10^-7 10^-8
$mu[0] = 0.000001; #min
$mu[1] = 0.0000001;
$mu[2] = 0.00000001; #max

# dump phylogeny only of clones that are above 10% of the neoplasm
$phylogeny_cutoff = 0.1;

# neighbor replacement probability
$p[0] = 0.0;
$p[1] = 1.0;
#$p[2] = 1.0;

#$p[1] = 0.25;
#$p[2] = 0.5;
#$p[3] = 0.75;
#$p[4] = 1.0;

# mutational model
$k[0] = 0;
$k[1] = 1;

# Count sweeps
$sweepcnt = 1;

for $ki (0..1) {
    for $pi (0..$#p){
	# Sweep selective mut rate
	for $mui (0..2){
	    # Sweep neutral mut rate
	    for $nui (0..0){
		for $rsmi (0..1){
		    # Repeat replicates times
		    for $replicates ($start_replicate..$stop_replicate){

		    open(out, ">infile$sweepcnt");
		    # K=1_P=0_MU=0.00001_NU=0.00001_RSM=1_B=0.02_D=0.016_NN=100_T=20_H=300_W=300
		    print out 
			"K=".$k[$ki]."_".
			"P=".$p[$pi]."_".
			"MU=".sprintf("%1.0e",$mu[$mui])."_".
			"NU=".sprintf("%1.0e",$nu[$nui])."_".
			"RSM=".$rsm[$rsmi]."_".
			"B=".$division_rate."_".
			"D=".$death_rate."_".
			"NN=".$neutral_loci_n."_".
			"T=".$simulation_t."_".
			"H=".$grid_size_y."_".
			"W=".$grid_size_x."\n";			
		    # Set the KnuthRandom seed, random int from 0 to 1 billion
		    print out int(rand(1000000000))."\n";
		    print out $replicates."\n"; # replicate number 1, 2, 3, ... so on
		    print out $grid_size_x."\n";
		    print out $grid_size_y."\n";
		    print out $neutral_loci_n."\n";

		    # get number of selective loci from RSM
		    @tokk = split(/\s/,$rsm_loci_string[$rsm[$rsmi]][1]);
		    print out ($#tokk+1)."\n";

		    # neighbor replacement prob
		    print out $p[$pi]."\n";
		    
		    # kappa
		    print out $k[$ki]."\n";

		    print out $division_rate."\n";
		    print out $division_max_rate."\n";
		    print out $division_min_rate."\n";

		    print out $death_rate."\n";
		    print out $death_max_rate."\n";
		    print out $death_min_rate."\n";

		    print out (sprintf("%.9f", $nu[$nui]))."\n";
		    print out (sprintf("%.9f", $nu[0]))."\n"; #fastest
		    print out (sprintf("%.9f", $nu[2]))."\n"; #slowest

		    print out (sprintf("%.9f", $mu[$mui]))."\n";
		    print out (sprintf("%.9f", $mu[0]))."\n"; #fastest
		    print out (sprintf("%.9f", $mu[2]))."\n"; #slowest
		    
		    print out $rsm_loci_string[$rsm[$rsmi]][1]."\n";
		    print out $rsm_loci_string[$rsm[$rsmi]][2]."\n";
		    print out $rsm_loci_string[$rsm[$rsmi]][3]."\n";

		    # Count time in generations (1) in days (0)		
		    print out "0\n";
		    print out $dump_grid_state_t."\n";
		    print out $cleanup_t."\n";
		    print out $simulation_t."\n";
		    # phylogeny
		    print out $phylogeny_cutoff."\n";
		    ## Seeding
		    print out q(250
0
0
300
0
0
0.1
0.0
0.0
);
		    
		    close(out);
		    		    			
			open(out2, ">job$sweepcnt.sh");
			print out2 "#!/bin/bash\n";		
			# Request 3.3 GB of memory per job
			#print out2 "java -server -Xmx2500m -cp /home/rumen/cryptsim/bin sim/DiversityTissueDriver infile$sweepcnt\n";
		    print out2 "java -server -Xmx2g -cp /home/rumen/cryptsim/bin sim/DiversityTissueDriver infile$sweepcnt\n";
			#print out2 "java -server -XX:+AggressiveHeap -cp /home/rumen/cryptsim/bin sim/DiversityTissueDriver infile$sweepcnt\n";
			close(out2);

		    # Increase count of sweeps
		    $sweepcnt++;

		    }

		}
	    }
	}
    }	
}
    

open(out,">submit.sge.pl");
print out q(#!/usr/bin/perl                                                                                                                                                        
	    if(defined($ARGV[1])){
		for($i=$ARGV[0]; $i<=$ARGV[1]; $i++){
		    $args = "qsub -l h_stack=8192K -q all.q -p -100 -cwd job$i.sh";
		    system($args);
		}
	    } else {
		$args = "qsub -l h_stack=8192K -q all.q -p -100 -cwd job".$ARGV[0].".sh";
		system($args);
	    }
);
close out;


